Logarithmic corrections in (4+l)-dimensional directed percolation 
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We simulate directed site percolation on two lattices with 4 spatial and 1 time-like dimensions 
(simple and body-centered hypercubic in space) with the standard single cluster spreading scheme. 
For efficiency, the code uses the same ingredients (hashing, histogram re-weighing, and improved 
estimators) as described in Phys. Rev. E 67, 036101 (2003). Apart from providing the most 
precise estimates for p c on these lattices, we provide a detailed comparison with the logarithmic 
corrections calculated by Janssen and Stenull [Phys. Rev. E 69, 016125 (2004)]. Fits with the 
leading logarithmic terms alone would give estimates of the powers of these logarithms which are 
too big by typically 50%. When the next-to-leading terms are included, each of the measured 
quantities (the average number of sites wetted at time t, their average distance from the seed, 
and the probability of cluster survival) can be fitted nearly perfectly. But these fits would not be 
mutually consistent. With a consistent set of fit parameters, one obtains still much improvement 
over the leading log - approximation. In particular we show that there is one combination of these 
three observables which seems completely free of logarithmic terms. 



I. INTRODUCTION 

Although it is well known that all critical phenom- 
ena have logarithmic corrections at their upper critical 
dimensions, and although the leading terms are easily 
calculated from the renormalization group, it is in gen- 
eral not easy to verify these predictions numerically. In 
equilibrium models, one reason is that it is difficult to 
simulate a sufficiently large system in high dimensions, 
both because of storage and of CPU requirements. The 
other reason is that, together with powers of the loga- 
rithm of the system size L, one usually has also terms 
of type log log L etc. If these are not known explicitly 
(and their computation is much more demanding), one 
has hardly any chance to verify the leading terms. 

The situation is somewhat better in models with long 
range interactions [lj and in tricritical phenomena (2J 
where the upper critical dimension is lower than in ordi- 
nary critical phenomena. It is also better in models like 
self-avoiding walks or percolation, where one does not 
need to simulate the entire lattice, but only fractal ob- 
jects with much lower dimension. For SAWs, e.g., it was 
possible to verify the structure of logarithmic corrections 
quite in detail [3[, since there one only has to simulate 
walks with dimension two, and since the next-to leading 
terms in the field theoretic treatment could be calculated. 

In the present paper we study directed percolation 
(DP). There, the upper critical dimension is 5. When in- 
terpreted as a spreading phenomenon, this corresponds to 
4 spatial dimensions. Critical clusters then have spatial 
fractal dimension Df — 2, i.e. then it becomes also feasi- 
ble to study systems with very large correlation lengths. 
In addition, the leading and next-to-leading logarithmic 
terms have been calculated recently from field theory Q , 
so that we have a good theoretical prediction to compare 



with. 

We study only site percolation, but on two lattices: 
the simple hypercubic (she) lattice in 4 dimension, and 
the body-centered hypercubic (bchc) lattice. The former 
has 2d = 8 neighbours which can be infected in each 
time step, the latter has 2 d = 16 neighbours. We use 
the standard spreading paradigm where we start with a 
single infected site and infect in each time step neigh- 
bouring sites with probability p. Sites stay infective for 
one time step, after that they become again susceptible. 
We measure the average number N(t) of infected sites, 
the r.m.s. distance R(t) of infected sites from the seed 
site, and the probability P(t) that there is still at least 
one infected site (i.e., that the cluster is still alive) at 
time t. The total sample sizes are 5.5 x 10 7 clusters for 
the she lattice, and 1.5 x 10 7 clusters for the bchc lattice, 
both with t max = 8000. 

The code used to simulate this is very similar to the 
one used in Q for high-dimensional ordinary percolation: 

1) We used hashing to store very large virtual lattices. 

2) In addition to the straightforward averages we also 
estimated in each run two averages obtained by re- 
weighing, corresponding to one p— value slightly above 
and to one p— value slightly below the point at which 
we simulate. This is equivalent to histogram re-weighing 
[f|0] j but avoids the need for storing huge histograms. 

3) We used improved estimators for N(t) and R(t), as 
described in [5j. These estimators were found to lead 
to large variance reduction (the same concept was used 
recently also for random walks with memory, where it 
also gave substantial improvements @). These estima- 
tors were found to lead to large variance reduction. Es- 
sentially, the idea is not to measure the actual number 
of offsprings in each generation (and their distances from 
the seed), but to measure the estimated number of off- 



2 



o 

CO 



-g 
c5 
■a 



V 

o 
o 

CD 



0.01 



0.001 



0.0001 



bcc, standard est. 

sc, standard est. 
bcc, improved est. 

sc, improved est. 



10 



100 
t 



1000 



10000 



FIG. 1: (color online) Log-log plot of statistical errors (one 
a) of N(t) against t. The upper curves are for the usual 
estimates, the lower ones are for the optimized improved es- 
timates. 
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FIG. 2: (color online) Average number of infected sites, N(t), 
for the she lattice. The three noisy curves are for p = p c and 
for p = p c ± Ap c . The other two curves show the leading log 
term (oc (ln^/to)) 1 ^ 6 with t' = 2) and the full prediction of 
Janssen and Stenull 3, Eq. (1) with to = 0.5, ti = 1.0. 



springs per active site (and their estimated distances). 
These estimates are made by counting the number of 
free neighbouring sites and multiplying it by p. This 
eliminates the fluctuations in the actual number of wet- 
ted sites resulting from the random number generator. 
Indeed, we found that the improved estimator gave not 
only smaller variances than the standard estimator, but 
that the covariances between the two happened to be 
negative (we have no explanation for this lucky coinci- 
dence). Thus we can optimize the estimator by taking 
that particular linear combination which has the small- 
est variance. The resulting errors for N(t) are shown in 
Fig.l. We see a reduction by roughly a factor 3, corre- 
sponding to a reduction of CPU time by a factor 10. The 
improvement was even larger (factor « 4) for R(t). For 
P(t) no similar improved estimator seems to exist. 



II. RESULTS 

Our main results are shown in Figs. 2 to 4. In each 
of them we show our results for the she lattice together 
with the leading logarithmic term and with a fit based on 
the full analytic results of The integration constants 
ti appearing in the logarithms are the same for all three 
observables. 

The results of [H can be rewritten as 
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x (1 + 0((lnlnt/lnt)^,lnlnt/ln t, 1/ln t) (1) 
with i = 1,2, and 3. Here, 

X 1 =N(i), X 2 =tP(t), X 3 = R 2 (t)/t, (2) 
the exponents cti are equal to 

a x = 1/6, a 2 = 1/2, a 3 = 1/12, (3) 
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FIG. 3: (color online) Cluster survival probability multiplied 
by t, tP(t), for the she lattice. The meaning of the curves is as 
for Fig. 2, except that the leading log term is oc (ln(t/i )) 1//2 
and that the three data curves are separated by 4Ap c . The 
values for io,£i, and t' are the same as in Fig. 2. 



the other known quantities are b = 1.30204, a x = 
0.1831, a 2 = -1.5193, a 3 = -1.7010, and t and t t 
are unknown integration constants from the renormaliza- 
tion group flow. Notice that to and t\ are not universal 
(they differ between models), but they are the same for 
all observables within one model - although using differ- 
ent values of ti for different observables could effectively 
take into account of higher order corrections. 

The first observation is that the leading logarithms 
alone are not sufficient to describe the data. Using only 
these terms, i.e. making ansatzes Xi = X^[la j^] ai , we 
would overestimate ap and ap by roughly 50%. The 
constant to can be chosen such that a nearly perfect fit 
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FIG. 4: (color online) Squared cluster radius divided by t, 
R 2 (t)/t, for the she lattice. The meaning of the curves is 
again as for Fig. 3, except that the leading log term is oc 
(ln(^)) 1/12 . 



0.04 
0.03 
0.02 
0.01 


-0.01 
-0.02 
-0.03 
-0.04 



bec, p = 0.07558474 

bec, p = 0.07558502 
bec, p = 0.07558530 
sc, p = 0.1461589 

sc, p = 0.1461592 
sc, p = 0.1461595 




10 



100 

t 



1000 10000 



FIG. 5: (color online) Difference between the average number 
of infected sites N(t) and optimized fits with Eq.(l), both for 
the she and for the bchc lattice. In the fit optimal values for 
to and ti are used, although these values would give poor fits 
when used for P(t) and R(t). 



is obtained for N(t) at large t. But this value of to gives 
bad results for the other two variables. Also, N(t) is the 
the variable which depends most sensitively on the exact 
value of p c . It is mainly for the latter that we need high 
statistics. Without a good estimate of p c we could not 
get a decent estimate of the logarithmic corrections from 
the leading terms alone. The same results were obtained 
for the bchc lattice (not shown here). Our estimates for 
p c are 

p c = 0.0755850 ± 0.0000003 (bchc), 

p c = 0.1461592 ±0.0000003 (she). (4) 

This is to be compared to Ref. [9], where the authors 
studied steady state DP with a weak rate h for "immi- 
gration" (i.e., sites are turned infective with a rate h, 
even when they have no infected neighbour), and then 
considered the limit h — > 0. The observable measured 
in [j| was the density of infected sites. Such simula- 
tions are of course much more cumbersome. In addition 
to corrections from the limit h — > one also has finite 
size corrections which are completely absent in spread- 
ing simulations. Indeed, the estimate for p c given in 0, 
p c = 0.075582 ± 0.000017 for the bchc lattice, has an er- 
ror about 60 times larger than ours. Nevertheless, very 
good agreement was found in Q when comparing with 
the leading log terms only. We believe that this is a bit 
fortuitous. 

The error estimates in Eq.(j4]) are of course subjective, 
as is true for all extrapolations and, in particular, also for 
any critical exponents. To support the above estimates 
we show in Fig. 5 our values of N(t), after having sub- 
tracting from them the best fits using Eq.{T]). In spite 
of the very small error bars of the raw data, the fits are 
perfect for t > 20. The lines seen in Fig. 5 correspond to 
p c ± Ap c , with Ape given in Eq.(|3]). 
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FIG. 6: (color online) Product P 2 (t)R A (t) / N(t) against t for 
the she lattice. Notice that the three curves correspond to 
values of p whose difference is about 30 standard deviations. 
The leading logarithmic correction would be linear in Int. 



Unfortunately the fits used in Fig. 5, although presum- 
ably correct for large values of t and therefore suitable 
for estimating p c , are not to be taken too seriously. This 
is seen from the fact that using the same values of to and 
t\ would give rather poor fits for the other two observ- 
ables. As a good compromise we used t n = 0.5, t\ = 1.0 
in Figs. 2 to 4. We see that none of the three fits is 
perfect, but all are quite reasonable and definitely give 
a big improvement over the leading term. Thus we can 
safely conclude that the field theoretic calculations of [3| 
are verified by our simulations. 

Equation (1) was indeed obtained in Ref. [4] by first 
deriving parametric forms X t = Xi(w) and t = t(w), 
and then inverting the latter to w — w(t). Since the 
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FIG. 7: (color online) Product Eq.(0 against t for the she 
lattice. The three data curves are for p c and for p c ± 4Ap c . 



parametric representations are only to lowest orders, the 
inversion introduces errors which, although subdominant 
asymptotically, might be numerically large. Comparing 
directly with the parametric expressions (Eqs.(15) and 
the first lines of Eqs.(25), (31), and (44) in Ref. [4]) gives 
indeed significant further improvements for small t. 

Before concluding, let us make two remarks. The first 
concerns hyperscaling. Usually, hyperscaling is formu- 
lated in terms of critical exponents. Writing N(t) ~ 
P(t) ~ t~ 5 , R 2 (t) ~ t z at p = p c , one expects for 
d < d c = 4 that dz/2 = 26 + rj. This is no longer true 
for d > 4 where rj = 0,8 = 1, and z = 2, but it still 
should hold in d = 4. Written in terms of the observ- 
ables themselves, hyperscaling is equivalent (for d < 4) 
to 

P 2 {t)R d {t) /N{t) w const. (5) 

From Eq. (JTJ) we sec that this should be violated by loga- 
rithmic terms at d = 4, 

P 2 {t)R d {t)/N{t) ~ [lnt]2«2+2a 3 -a 1 = (g) 

We see from Fig. 6 that this product indeed increases 
strongly with t, but the increase is far from linear in hit. 
Thus, next-to-leading terms again are important. The 
corrections given in Eq.(l) give a big improvement, al- 
though they are not perfect. An interesting observation 



is that this product depends very weakly on p, making 
it thus an ideal test object for further non-leading loga- 
rithmic corrections. 

The second remark concerns another product of 
N(t),P(t), and R(t). Using Eq.([T]) we can form one 
combination (and of course all its powers) which con- 
tains, up to the order considered in Eq.|T]), no loga- 
rithmic corrections at all. It is given by ]X with 
^Zif'iOii = 'Yliit i 'i 0i i a i — 0- Numerically, we thus obtain 
that 




10.8427 

/ (tP(t)) « const. (7) 



We plot this combination in Fig. 7, together with a fit of 
the type a + b/t A . Numerically we found A = 0.75. 
Of course one should not take this exponent very se- 
rious (it could well be that the correct exponent is 
1/2 or 1), but it seems rather convincing that logarith- 
mic terms are completely absent. Notice that this is 
not trivial. A priori, we should have expected terms 
~ 0((lnlnt/lnt) 2 ,lnlnt/ln 2 t,l/ln 2 t). This might hint 
at a special structure of the renormalization group flow, 
although this does not seem likely from the way in which 
Eq. (1) was derived [1(3]. 



III. SUMMARY 

We have shown that improved algorithms for cluster 
spreading allow, even with rather modest effort (the to- 
tal CPU time used for this paper was about 1 week on 
a fast PC), a rather stringent verification of logarithmic 
corrections at the upper critical dimension of one of the 
standard non-equilibrium critical phenomena. A prereq- 
uisite for this was, however, the availability of more than 
the leading log terms. If we would have had only the 
leading terms available for comparison (as was the case 
for the steady-state equation of state studied in 0), even 
with much more CPU time only an order of magnitude 
verification would have been possible. 
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